﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Numerics;
using MathNet.Numerics;//需调用开源库MathNet.Numerics

using CC.Extending;
using System.Diagnostics;
using System.Text.RegularExpressions;
using System.Runtime.CompilerServices;
using CC.Common;

namespace CcMath
{
    public enum ECoordinates
    {
        Cartesian,Cyclinder,Spherical
    }
    
    
    /// <summary>
    /// 引用了MathNet.Numerics 开源数学库，功能更强大
    /// </summary>
    public class MyMath
    {
        private MyMath() { }
       
        /// <summary>
        /// 求解二维数据矩阵的差值数组：dataArry中，每组数据的最大差值
        /// 特别用于求解端口差异值
        /// 2020-10-4
        /// </summary>
        /// <param name="dataArray">二维数组，要求每个数组的长度相同。工程中，往往是行是端口，列是频率</param>
        /// <param name="correctPeriod">是否左周期性矫正。由于相位存在周期性，需要做矫正</param>
        /// <returns></returns>
        public static double[] CalculateRangeArray(IEnumerable<double[]> dataArray, bool correctPeriod)
        {
            //检查列表中每列的数据长度是否相同，如果不同，则报错
            int length = dataArray.FirstOrDefault().Count();
            foreach (var item in dataArray)
            {
                if (item.Count() != length)
                {
                    throw (new Exception("输入的数据长度不一致，无法做数据差值运算。\r\n可能的原因是起止频率或频率点数不一致，请检查数据"));
                }
            }
            //计算S参数每个频率的端口差异
            double[][] transSP = dataArray.ToArray().Transpose();//转置,变成每个频率的不同端口值（行为频率，列为端口）
            double[] ranges = new double[transSP.Length];
            //求解范围矩阵
            for (int f = 0; f < ranges.Length; f++)
            {
                var delta = transSP[f].Max() - transSP[f].Min();
                if (correctPeriod)
                {
                    //考虑相位的周期性跳变，对delta进行矫正：取delta，delta+/-360的最小绝对值
                    var datas = new double[] { Math.Abs(delta), Math.Abs(delta + 360f), Math.Abs(delta - 360f) };
                    ranges[f] = datas.Min();
                }
                else
                {
                    ranges[f] = delta;
                }
            }
            return ranges;
        }
        /// <summary>
        /// 数值积分
        /// </summary>
        /// <param name="pointsAscend">按x排序</param>
        /// <param name="startIndex">起始索引，默认为0，表示从头开始积分</param>
        /// <param name="endIndex">终点索引，默认为0，表示计算到末尾</param>
        /// <returns></returns>
        public static double Integrate(Point32[]pointsAscend,int startIndex=0,int endIndex=0)
        {
            double res = 0;
            if (endIndex == 0) endIndex = pointsAscend.Length - 1;
            if(startIndex<endIndex)
            {
                for (int i = startIndex; i < endIndex; i++)
                {
                    double deltaX = pointsAscend[i + 1].X - pointsAscend[i].X;
                    double y = (pointsAscend[i].Y + pointsAscend[i + 1].Y) / 2;
                    double num = deltaX * y;
                    res += num;
                }
            }
            
            return res;
        }
        /// <summary>
        /// 寻找数列中的所有0值点
        /// </summary>
        /// <param name="points"></param>
        /// <param name="order">阶数，意思是取左右各多少个点来做判断</param>
        /// <returns>0值点的索引清单</returns>
        public static List<int> FindZeroIndex(Point32[] points,int order=2)
        {
            
            List<int> list = new List<int>();
            for (int i = 0; i < points.Length; i++)
            {
                var subarr = points.SubArrayPeriod(i - order, i + order + 1);
                bool hasZero = true;//预设为是0点；x[i]*x[length-i]<0,表示两个之间有0值点
                for (int j = 0; j < order; j++)
                {
                    bool res = subarr[j].Y * subarr[subarr.Length - j - 1].Y < 0;
                    hasZero &= res;
                }
                
                if (hasZero)
                {
                    list.Add(i);
                }
            }
            return list;
        }
        /// <summary>
        /// 搜索0点（极小值点）
        /// </summary>
        /// <param name="points"></param>
        /// <param name="order">阶数，意思是取左右各多少个点来做判断</param>
        /// <returns>极小值点清单</returns>
        public static List< Point32>FindNulls(Point32[] points, int order = 2)
        {
            List<Point32> list = new List<Point32>();
            for (int i = 0; i < points.Length; i++)
            {
                var subarr = points.SubArrayPeriod(i - order, i + order + 1);
                var res = subarr.Where(a => a.X != points[i].X && a.Y <= points[i].Y);
                if (res.Count() == 0)//左右两边的y均小于中心值
                {
                    list.Add(points[i]);
                }
            }

            return list;
        }

        /// <summary>
        /// 计算球坐标下给定方向矢量在给定参考方向上投影长度
        /// phi从x轴正向开始，sita从z轴正向开始
        /// </summary>
        /// <param name="vector">方向矢量[r, phi, sita]</param>
        /// <param name="refer">参考方向[r, phi, sita]</param>
        /// <returns></returns>
        public static double SphericalShadow(Point3D vector, Point3D refer)
        {
            double cosPsai = Math.Sin(refer.Sita.ToRad()) * Math.Sin(vector.Sita.ToRad()) * Math.Cos((vector.PhiofSph - refer.PhiofSph).ToRad())
                + Math.Cos(vector.Sita.ToRad()) * Math.Cos(refer.Sita.ToRad());
            var shadow = vector.Magtitude * cosPsai;
            return shadow;
        }
        /// <summary>
        /// 将输入的值限制在指定的周期内
        /// </summary>
        /// <param name="origin">原始值</param>
        /// <param name="periodicity">周期数</param>
        /// <param name="minValue">周期内最小值</param>
        /// <returns></returns>
        public static float LimitToPeriod(float origin, float periodicity, float minValue = 0)
        {
            if (origin < minValue)
            {
                origin += periodicity;
                return LimitToPeriod(origin, periodicity, minValue);
            }
            else if (origin > minValue + periodicity)
            {
                origin -= periodicity;
                return LimitToPeriod(origin, periodicity, minValue);
            }
            else
            {
                return origin;
            }

        }
        #region 数据平滑算法
        /// <summary>
        /// 使用七点平滑算法对数据点进行平滑处理
        /// </summary>
        /// <param name="points"></param>
        /// <param name="num"></param>
        public static Point32[] SmoothPoints7(Point32[] points, int num)
        {
            float[] data = new float[points.Length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = points[i].Y;
            }
            Smooth(data, num);
            for (int i = 0; i < points.Length; i++)
            {
                points[i].Y = data[i];
            }
            return points;
        }
        /// <summary>
        /// 使用七点平滑算法对数据点进行平滑处理
        /// </summary>
        /// <param name="points"></param>
        /// <param name="num"></param>
        public static List<Point32> SmoothPoints7(List< Point32> points, int num)
        {
            float[] data = new float[points.Count];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = points[i].Y;
            }
            Smooth(data, num);
            Point32[] dps = new Point32[data.Length];
            for (int i = 0; i < points.Count; i++)
            {
                dps[i].Y = data[i];
                dps[i].X = points[i].X;
            }
            points = dps.ToList();
            return points;
        }
        /// <summary>
        /// 通过七点平滑算法对数据数组进行平滑处理
        /// </summary>
        /// <param name="x">要平滑的数据</param>
        /// <param name="num">平滑次数,5到10次效果比较明显</param>
        public static double[] Smooth7(double[] x, int num)
        {
            int length = x.Length;
            if (length >= 7 && num > 0)
            {
                for (int p = 1; p <= num; p++)
                {
                    for (int q = 0; q < x.Length; q++)
                    {
                        x[q] = (1 / 42.0) * (-4 * x.Fetch(q - 2) + 6 * x.Fetch(q - 1) + 12 * x.Fetch(q) + 14 * x.Fetch(q + 1) + 12 * x.Fetch(q + 2) + 6 * x.Fetch(q + 3) - 4 * x.Fetch(q + 4));
                    }
                    //修改阵列单元提取算法，索引不再受限 2021-4-23
                    //x[0] =  (31 * x[0] + 8 * x[1] - 4 * x[2] - 4 * x[3] + x[4]+4*x[5]-2*x[6])/42f;
                    //x[1] =  (8 * x[0] + 19 * x[1] +16 * x[2] +6 * x[3] -4* x[4] -7 * x[5] +4 * x[6])/42f;
                    //x[2] = (-4 * x[0] + 16 * x[1] + 19 * x[2] + 12 * x[3] +2 * x[4] - 4 * x[5] +  x[6])/42f;
                    //for (int q = 3; q < length - 3; q++)
                    //{
                    //    x[q] = (1 / 42.0) * (-4 * x[q - 2] + 6 * x[q - 1] + 12 * x[q] + 14 * x[q + 1] +12 * x[q + 2]+6*x[q+3]-4*x[q+4]);
                    //}
                    //x[length - 3] = ( x[length - 7] -4 * x[length - 6] + 2 * x[length - 5] + 12 * x[length - 4] + 19 * x[length - 3]+16*x[length-2]-4*x[length-1])/42f;
                    //x[length - 2] = (4*x[length - 7] - 7 * x[length - 6] -4 * x[length - 5] + 6 * x[length - 4] + 16 * x[length - 3] + 19 * x[length - 2]+ 8 * x[length - 1])/42f;
                    //x[length - 1] = (-2 * x[length - 7] +4 * x[length - 6]+ x[length - 5] -4 * x[length - 4] -4 * x[length - 3] + 8 * x[length - 2] +39 * x[length - 1]) / 42f;
                }
            }
            return x;
        }

       

        /// <summary>
        /// 通过七点平滑算法对数据数组进行平滑处理
        /// </summary>
        /// <param name="x">要平滑的数据</param>
        /// <param name="num">平滑次数,5到10次效果比较明显</param>
        public static void Smooth7(float[] x, int num)
        {
            int length = x.Length;
            if (length >= 7 && num > 0)
            {
                for (int p = 1; p <= num; p++)
                {
                    for (int q = 0; q < x.Length; q++)
                    {
                        x[q] = (float)(1 / 42.0) * (-4 * x.Fetch(q - 2) + 6 * x.Fetch(q - 1) + 12 * x.Fetch(q) + 14 * x.Fetch(q + 1) + 12 * x.Fetch(q + 2) + 6 * x.Fetch(q + 3) - 4 * x.Fetch(q + 4));
                    }
                    //修改阵列单元提取算法，索引不再受限 2021-4-23
                    //x[0] = (31 * x[0] + 8 * x[1] - 4 * x[2] - 4 * x[3] + x[4] + 4 * x[5] - 2 * x[6]) / 42f;
                    //x[1] = (8 * x[0] + 19 * x[1] + 16 * x[2] + 6 * x[3] - 4 * x[4] - 7 * x[5] + 4 * x[6]) / 42f;
                    //x[2] = (-4 * x[0] + 16 * x[1] + 19 * x[2] + 12 * x[3] + 2 * x[4] - 4 * x[5] + x[6]) / 42f;
                    //for (int q = 3; q < length - 3; q++)
                    //{
                    //    x[q] = (float)(1 / 42.0) * (-4 * x[q - 2] + 6 * x[q - 1] + 12 * x[q] + 14 * x[q + 1] + 12 * x[q + 2] + 6 * x[q + 3] - 4 * x[q + 4]);
                    //}
                    //x[length - 3] = (x[length - 7] - 4 * x[length - 6] + 2 * x[length - 5] + 12 * x[length - 4] + 19 * x[length - 3] + 16 * x[length - 2] - 4 * x[length - 1]) / 42f;
                    //x[length - 2] = (4 * x[length - 7] - 7 * x[length - 6] - 4 * x[length - 5] + 6 * x[length - 4] + 16 * x[length - 3] + 19 * x[length - 2] + 8 * x[length - 1]) / 42f;
                    //x[length - 1] = (-2 * x[length - 7] + 4 * x[length - 6] + x[length - 5] - 4 * x[length - 4] - 4 * x[length - 3] + 8 * x[length - 2] + 39 * x[length - 1]) / 42f;
                }
            }
        }
        /// <summary>
        /// 使用五点二次平滑算法对数据点进行平滑处理
        /// </summary>
        /// <param name="points"></param>
        /// <param name="num"></param>
        public static Point32[] SmoothPoints(Point32[] points, int num)
        {
            float[] data = new float[points.Length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = points[i].Y;
            }
            data = Smooth(data, num);
            Point32[] pts = new Point32[points.Length];
            for (int i = 0; i < pts.Length; i++)
            {
                pts[i] = new Point32(points[i].X, data[i]);
            }
            return pts;
        }
        /// <summary>
        /// 通过五点二次平滑算法对数据数组进行平滑处理
        /// </summary>
        /// <param name="dataArray">要平滑的数据</param>
        /// <param name="num">平滑次数,5到10次效果比较明显</param>
        public static float[] Smooth(float[] dataArray, int num)
        {
            int length = dataArray.Length;
            if (dataArray.Length >= 5 && num > 0)
            {
                for (int p = 1; p <= num; p++)
                {
                   
                    for (int q = 0; q < dataArray.Length; q++)
                    {
                        dataArray[q]= Convert.ToSingle(1 / 35.0) * (-3 * dataArray.Fetch(q - 2) + 12 * dataArray.Fetch(q - 1) + 17 * dataArray.Fetch(q) + 12 * dataArray.Fetch(q + 1) - 3 * dataArray.Fetch(q + 2));
                    }
                    //修改阵列单元提取算法，索引不再受限 2021-4-23
                    //dataArray[0] = Convert.ToSingle(1 / 35.0) * (31 * dataArray[0] + 9 * dataArray[1] - 3 * dataArray[2] - 5 * dataArray[3] + 3 * dataArray[4]);
                    //dataArray[1] = Convert.ToSingle(1 / 35.0) * (9 * dataArray[0] + 13 * dataArray[1] + 12 * dataArray[2] + 6 * dataArray[3] - 5 * dataArray[4]);
                    //for (int q = 2; q < dataArray.Length - 2; q++)
                    //{
                    //    dataArray[q] = Convert.ToSingle(1 / 35.0) * (-3 * dataArray[q - 2] + 12 * dataArray[q - 1] + 17 * dataArray[q] + 12 * dataArray[q + 1] - 3 * dataArray[q + 2]);
                    //}
                    ////从SLD处拿到的代码没有下面这两行
                    //dataArray[length - 2] = (1 / 35f) * (-5 * dataArray[length - 5] + 6 * dataArray[length - 4] + 12 * dataArray[length - 3] + 13 * dataArray[length - 2] + 9 * dataArray[length - 1]);
                    //dataArray[length - 1] = (1 / 35f) * (3 * dataArray[length - 5] + 5 * dataArray[length - 4] - 3 * dataArray[length - 3] + 9 * dataArray[length - 2] + 31 * dataArray[length - 1]);
                }

            }
            return dataArray;
        }

        /// <summary>
        /// 通过五点二次平滑算法对数据数组进行平滑处理
        /// </summary>
        /// <param name="dataArray">要平滑的数据</param>
        /// <param name="num">平滑次数,5到10次效果比较明显</param>
        public static double[] Smooth(double[] dataArray, int num)
        {
            int length = dataArray.Length;
            if (length >= 5 && num > 0)
            {
                for (int p = 1; p <= num; p++)
                {
                    for (int q = 0; q < dataArray.Length; q++)
                    {
                        dataArray[q] = (1 / 35.0) * (-3 * dataArray.Fetch(q - 2) + 12 * dataArray.Fetch(q - 1) + 17 * dataArray.Fetch(q) + 12 * dataArray.Fetch(q + 1) - 3 * dataArray.Fetch(q + 2));
                    }
                    //修改阵列单元提取算法，索引不再受限 2021-4-23
                    //dataArray[0] = (1 / 35.0) * (31 * dataArray[0] + 9 * dataArray[1] - 3 * dataArray[2] - 5 * dataArray[3] + 3 * dataArray[4]);
                    //dataArray[1] = (1 / 35.0) * (9 * dataArray[0] + 13 * dataArray[1] + 12 * dataArray[2] + 6 * dataArray[3] - 5 * dataArray[4]);
                    //for (int q = 2; q < length - 2; q++)
                    //{
                    //    dataArray[q] = (1 / 35.0) * (-3 * dataArray[q - 2] + 12 * dataArray[q - 1] + 17 * dataArray[q] + 12 * dataArray[q + 1] - 3 * dataArray[q + 2]);
                    //}
                    ////从SLD处拿到的代码没有下面这两行
                    //dataArray[length - 2] = (1 / 35f) * (-5 * dataArray[length - 5] + 6 * dataArray[length - 4] + 12 * dataArray[length - 3] + 13 * dataArray[length - 2] + 9 * dataArray[length - 1]);
                    //dataArray[length - 1] = (1 / 35f) * (3 * dataArray[length - 5] + 5 * dataArray[length - 4] - 3 * dataArray[length - 3] + 9 * dataArray[length - 2] + 31 * dataArray[length - 1]);
                }
            }
            return dataArray;
        }

        /// <summary>
        /// 通过五点二次平滑算法对数据数组进行平滑处理
        /// </summary>
        /// <param name="dataArray">要平滑的数据</param>
        /// <param name="num">平滑次数,5到10次效果比较明显</param>
        public static Complex[] Smooth(Complex[] dataArray, int num)
        {
            double[] real = new double[dataArray.Length];
            double[] imag = new double[dataArray.Length];
            for (int i = 0; i < dataArray.Length; i++)
            {
                real[i] = dataArray[i].Real;
                imag[i] = dataArray[i].Imaginary;
            }
            Smooth(real, num);
            Smooth(imag, num);
            for (int i = 0; i < dataArray.Length; i++)
            {
                dataArray[i] = new Complex(real[i], imag[i]);
            }
            return dataArray;
        }
        #endregion

        /// <summary>
        /// 用一元二次方程拟合数据点，并求解其根
        /// </summary>
        /// <param name="points">三个点，跨越0点的三个点。不足抛异常，多于则取前三</param>
        /// <returns>在points范围内的拟合方程的根</returns>
        public static double FindZeroWithQuadraticFit(Point32[] points)
        {
            if(points.Length<3) throw (new Exception("一元二次方程拟合，至少必须三个数据点！"));
            (double a, double b, double c) = EvalCoefficientOfQuadraticEquation(points);
            (double x1, double x2) = EvalRootOfQuadraticEquation(a, b, c);
            if (x1 >= points.FirstOrDefault().X && x1 <= points.Last().X) return x1;
            else if (x2 >= points.FirstOrDefault().X && x2 <= points.Last().X) return x2;
            else 
            { 
                throw (new Exception("找不到points区间内的0点！")); 
            }
        }
        /// <summary>
        /// 求解一元二次方程的系数
        /// </summary>
        /// <param name="points">三个点，不足抛异常，多于则取前三。</param>
        /// <returns></returns>
        private static (double a,double b,double c) EvalCoefficientOfQuadraticEquation(Point32[] points)
        {
            if(points.Length<3)
            {
                throw (new Exception("求一元二次方程系数：需要三个点，但仅传入了" + points.Length + "个"));
            }
            double fenzi = (points[2].Y - points[1].Y) * (points[1].X - points[0].X) - (points[1].Y - points[0].Y) * (points[2].X - points[1].X);
            double fenmu1 = (points[2].X * points[2].X - points[1].X * points[1].X) * (points[1].X - points[0].X);
            double fenmu2= (points[1].X * points[1].X - points[0].X * points[0].X) * (points[2].X - points[1].X);
            double a = fenzi / (fenmu1 - fenmu2);
            fenzi = (points[1].Y - points[0].Y) - a * (points[1].X * points[1].X - points[0].X * points[0].X);
            double b = fenzi / (points[1].X - points[0].X);
            double c = points[0].Y - a * points[0].X * points[0].X - b * points[0].X;
            return (a, b, c);
        }
        /// <summary>
        /// 求解一元二次方程y = ax^2 + bx + c的两个根
        /// </summary>
        /// <param name="a">y = ax^2 + bx + c的a</param>
        /// <param name="b">y = ax^2 + bx + c的b</param>
        /// <param name="c">y = ax^2 + bx + c的c</param>
        /// <returns>(x1-负根，x2-正根)</returns>
        private static (double x1, double x2) EvalRootOfQuadraticEquation(double a,double b,double c)
        {
            double x1, x2;
            if(a==0)//退化成一元一次方程
            {
                x1 = -c / b;
                x2 = x1;
            }
            else
            {
                x1 = (-b + Math.Sqrt(b * b - 4 * a * c)) / (2 * a);
                x2 = (-b - Math.Sqrt(b * b - 4 * a * c)) / (2 * a);
            }

            
            return (x1, x2);
        }
        /// <summary>
        /// 一维插值产生新的数据点
        /// </summary>
        /// <param name="srcPoints">源数据点,需按X轴数据从小到大排列</param>
        /// <param name="x">数据点自变量</param>
        /// <returns></returns>
        public static Point32 Interpolate1(Point32[] srcPoints,double x)
        {
            
            double[] xArray = new double[srcPoints.Length];
            double[] yArray = new double[srcPoints.Length];
            for (int i = 0; i < xArray.Length; i++)
            {
                xArray[i] = srcPoints[i].X;
                yArray[i] = srcPoints[i].Y;
            }
            double y = Interpolate1(yArray, xArray, x);
            return new Point32((float)x, (float)y);
        }
        /// <summary>
        /// 一维插值函数
        /// 已知xArray,yArray,求指定点x处的y值
        /// </summary>
        /// <param name="yArray">y向量</param>
        /// <param name="xArray">对应的x坐标向量，从小到大排列</param>
        /// <param name="x">指定x值</param>
        /// <returns></returns>
        public static double Interpolate1(float[] yArray, float[] xArray, double x)
        {
            if (yArray.Length != xArray.Length)
            {
                throw (new Exception("xArray和yArray的长度必须相等！"));
            }
            if (yArray.Length < 2)
            {
                throw (new Exception("xArray和yArray的长度不可以小于2！"));
            }
            //如果x包含在xArray中，则不必插值，直接返回值
            int xIndex = Array.IndexOf(xArray, x);
            if (xIndex >= 0) return yArray[xIndex];

            int dataLength = yArray.Length;
            double x1 = xArray[0], x2 = xArray[dataLength - 1], y1 = yArray[0], y2 = yArray[dataLength - 1];
            if (x < xArray[0])
            {
                x1 = xArray[0]; x2 = xArray[1];
                y1 = yArray[0]; y2 = yArray[1];
            }
            else if (x > xArray[dataLength - 1])
            {
                x1 = xArray[dataLength - 2]; x2 = xArray.Last();
                y1 = yArray[dataLength - 2]; y2 = yArray.Last();
            }
            else
            {

                for (int i = 0; i < dataLength; i++)
                {
                    if (xArray[i] > x)
                    {
                        x1 = xArray[i - 1]; x2 = xArray[i];
                        y1 = yArray[i - 1]; y2 = yArray[i];
                        break;
                    }
                }
            }

            return Interpolate1(x, x1, x2, y1, y2);
        }
        /// <summary>
        /// 一维插值函数
        /// 已知xArray,yArray,求指定点x处的y值
        /// </summary>
        /// <param name="yArray">y向量</param>
        /// <param name="xArray">对应的x坐标向量，从小到大排列</param>
        /// <param name="x">指定x值</param>
        /// <returns></returns>
        public static double Interpolate1(double[] yArray, double[] xArray, double x)
        {
            if (yArray.Length != xArray.Length)
            {
                throw (new Exception("xArray和yArray的长度必须相等！"));
            }
            if (yArray.Length < 2)
            {
                throw (new Exception("xArray和yArray的长度不可以小于2！"));
            }
            //如果x包含在xArray中，则不必插值，直接返回值
            int xIndex = Array.IndexOf(xArray, x);
            if (xIndex >= 0) return yArray[xIndex];

            int dataLength = yArray.Length;
            double x1 = xArray[0], x2 = xArray[dataLength - 1], y1 = yArray[0], y2 = yArray[dataLength - 1];
            if (x < xArray[0])
            {
                x1 = xArray[0]; x2 = xArray[1];
                y1 = yArray[0]; y2 = yArray[1];
            }
            else if (x > xArray[dataLength - 1])
            {
                x1 = xArray[dataLength-2]; x2 = xArray.Last();
                y1 = yArray[dataLength-2]; y2 = yArray.Last();
            }
            else
            {
               
                for (int i = 0; i < dataLength; i++)
                {
                    if ( xArray[i]>x )
                    {
                        x1 = xArray[i-1]; x2 = xArray[i];
                        y1 = yArray[i-1]; y2 = yArray[i];
                        break;
                    }
                }
            }

            return Interpolate1(x, x1, x2, y1, y2);
        }
        public static double Interpolate1(double x,double x1,double x2,double y1,double y2)
        {
           
            double y = y1 + (y2 - y1) * (x - x1) / (x2 - x1);
            return y;
        }

        /// <summary>
        /// 复数一维线性插值程序
        /// </summary>
        /// <param name="y"></param>
        /// <param name="x"></param>
        /// <param name="x1"></param>
        /// <returns></returns>
        public static Complex Interpolate1(Complex[] y, double[] x, double xn)
        {
            int index = -1;
            double x1 = 0, x2 = 0;
            Complex y1 = Complex.Zero, y2 = Complex.Zero, yn = Complex.Zero;
            if (xn < x[0]) //插值点在参考坐标左边,变成外推法
            {
                x1 = x[0]; x2 = x[1];
                y1 = y[0]; y2 = y[1];
            }
            else if (xn > x[x.Length - 1])//插值点在参考坐标右边，变成外推法
            {
                x1 = x[x.Length - 1]; x2 = x[x.Length - 2];
                y1 = y[y.Length - 1]; y2 = y[y.Length - 2];
            }
            else
            {
                for (int n = 0; n < x.Length; n++)
                {
                    if (x[n] == xn) index = n;
                    if (n != 0)
                    {
                        if (x[n - 1] < xn && x[n] > xn)
                        {
                            x1 = x[n - 1]; x2 = x[n];
                            y1 = y[n - 1]; y2 = y[n];
                        }
                    }
                }
            }
            if (index != -1) yn = y[index];//数值已含，不需插值
            else
            {
                yn = Interpolate1(xn, x1, x2, y1, y2);
            }
            return yn;
        }
        public static Complex Interpolate1(double x, double x1, double x2, Complex y1, Complex y2)
        {
            double k = (y1.Real - y2.Real) / (x2 - x1);
            double realPart = y1.Real + k * (x - x1);
            k = (y1.Imaginary - y2.Imaginary) / (x2 - x1);
            double imagPart = y1.Imaginary + k * (x - x1);
            return new Complex(realPart, imagPart);
        }
        public static double[] Interpolate1(double[,] yArray, double[] xAxis, double x)
        {
            int cols = yArray.GetLength(1);
            double[] iValue = new double[cols];
            for (int i = 0; i < cols; i++)
            {
                int rows = yArray.GetLength(0);
                double[] yArr = new double[rows];
                for (int r = 0; r < rows; r++)
                {
                    yArr[r] = yArray[r, i];
                }
                iValue[i] = Interpolate1(yArr, xAxis, x);
            }
            return iValue;
        }
        /// <summary>
        /// 求数组中最接近给定值的索引
        /// </summary>
        /// <param name="array">数组</param>
        /// <param name="value">给定值</param>
        /// <returns>最接近给定值的索引</returns>
        public static int IndexOf(Point32[] array, double xValue)
        {
            double[] ar = new double[array.Length];
            for (int i = 0; i < ar.Length; i++)
            {
                ar[i] = Math.Abs(array[i].X - xValue);
            }
            double min = ar.Min();
            for (int i = 0; i < ar.Length; i++)
            {
                ar[i] -= min;
            }
            for (int i = 0; i < ar.Length; i++)
            {
                if (ar[i] == 0) return i;
            }
            return -1;
        }
        /// <summary>
        /// 求数组中最接近给定值的索引
        /// </summary>
        /// <param name="array">数组</param>
        /// <param name="value">给定值</param>
        /// <returns>最接近给定值的索引</returns>
        public static int IndexOf(double[] array, double value)
        {
            double[] ar = new double[array.Length];
            for (int i = 0; i < ar.Length; i++)
            {
                ar[i] = Math.Abs(array[i] - value);
            }
            double min = ar.Min();
            for (int i = 0; i < ar.Length; i++)
            {
                ar[i] -= min;
            }
            for (int i = 0; i < ar.Length; i++)
            {
                if (ar[i] == 0) return i;
            }
            return -1;
        }
        
        /// <summary>
        /// 寻找峰值
        /// </summary>
        /// <param name="array"></param>
        /// <returns>峰值清单，如果没有峰值则返回空值</returns>
        /// <param name="waveletFilter">要过滤到的纹波</param>
        public static double[] FindPeaks(double[] array, float waveletFilter = 0)
        {
            if (array.Length < 1)
            {
                throw (new Exception("Array length Must bigger than 1"));
            }
            List<double> peaksL = new List<double>();

            for (int i = 1; i < array.Length - 1; i++)
            {
                if (array[i] >= array[i - 1]+ waveletFilter && array[i] >= array[i + 1]+ waveletFilter)
                {
                    peaksL.Add(array[i]);
                }
            }
            if (peaksL.Count <= 0)
            {
                return null;                
            }            
            return peaksL.ToArray();
        }
        /// <summary>
        /// 查找所有峰值。
        /// </summary>
        /// <param name="array"></param>
        /// <param name="peakIndexs">峰值点清单，如果没有则空值</param>
        /// <returns>峰值清单，如果没有峰值则返回第一个值</returns>
        public static double[] FindPeaks(double[] array, out int[] peakIndexs ,float waveletFilter = 0)
        {
            if (array.Length < 1)
            {
                throw (new Exception("Array length Must bigger than 1"));
            }
            List<double> peaksL = new List<double>();
            List<int> indList = new List<int>();
            for (int i = 1; i < array.Length - 1; i++)
            {
                if (array[i] >= array[i - 1] + waveletFilter && array[i] >= array[i + 1] + waveletFilter)
                {
                    peaksL.Add(array[i]);
                    indList.Add(i);
                }
            }
            if(peaksL.Count<=0)
            {//如果找不到峰值，就输出空值
                peakIndexs = null;
                return null;
            }
            else
            {
                peakIndexs = indList.ToArray();
            }
            return peaksL.ToArray();
        }
        /// <summary>
        /// 用多项式拟合方法找指定区间内阵列的最大值
        /// 注1：找出来的最大值是拟合出来的，通常与数列的最大值不相同
        /// 注2：提供的曲线最好是单峰的，元素数量3-5个为宜，否则拟合效果会比较差。
        /// </summary>
        /// <param name="yArray">数组</param>
        /// <param name="xArray">对应yArray的自变量数组</param>
        /// <param name="x">输出最大值自变量</param>
        /// 
        /// <returns>曲线拟合和的最大值</returns>
        public static double FindMaxCurveFit(double[] yArray, double[] xArray, out double xMax, double xStep = 0.01)
        {
            if(xArray.Length<= yArray.Length)
            {
                throw (new Exception("多项式拟合法求数组最大值数据异常：x数组和y数组的长度不匹配!"));
            }
            var yMax = yArray.Max();
            int indexMax = Array.IndexOf(yArray, yMax);
            xMax = xArray[indexMax];
            if (indexMax == 0 || indexMax == yArray.Length - 1) return yMax;
            double[] x;
            double[] y;
            if (indexMax >= 2 && indexMax < yArray.Length - 2)
            {
                y = new double[] { yArray[indexMax - 2], yArray[indexMax - 1], yArray[indexMax], yArray[indexMax + 1], yArray[indexMax + 2] };
                x = new double[] { xArray[indexMax - 2], xArray[indexMax - 1], xArray[indexMax], xArray[indexMax + 1], xArray[indexMax + 2] };
            }
            else
            {
                y = new double[] {  yArray[indexMax - 1], yArray[indexMax], yArray[indexMax + 1] };
                x = new double[] {  xArray[indexMax - 1], xArray[indexMax], xArray[indexMax + 1] };
            }
            double[] p = Fit.Polynomial(x,y, y.Length - 1);
            double deltaX = (x.Max() - x.Min()) / (100f);//步长缩短100倍去搜索最大值
            deltaX = deltaX < xStep ? deltaX : xStep;
            int N = (int)((x.Max() - x.Min()) / deltaX);
            double[] yFitted = new double[N];
            for (int i = 0; i < N; i++)
            {
                xMax = x.Min() + i * deltaX;
                yFitted[i] = PolyValue(p, xMax);
            }
            double ymax = yFitted.Max();
            return ymax;
        }
        /// <summary>
        /// 用多项式拟合方法找指定位置附近的阵列最大点
        /// 注1：找出来的最大值是拟合出来的，通常与数列的最大值不相同
        /// 注2：以-1dBc判断是否栅瓣
        /// </summary>
        /// <param name="points">要搜索的数组</param>
        /// <param name="nearby">参考点</param>
        /// <param name="xStep">x搜索步长</param>
        /// <returns></returns>
        public static Point32 FindMaxNearTo_CurveFit(Point32[] points, double nearby, double xStep = 0.01)
        {
              var count = points.Count(pp => float.IsNaN(pp.Y));
            if (count > 0)
            {
                throw new Exception("输入数据含有NaN值，无法参与计算！");
            }
            //查找看看是否有栅瓣：有多个大于-1dBc的峰值点

            var peaks = MyMath.FindPeaks(points);
            var gratings = peaks.Where(d => d.Y >= peaks.Max(p => p.Y) - 1).ToArray();
            if (gratings.Count() > 1)//有栅瓣，取最接近nearby的那个波瓣作为主瓣
            {
                var nears = gratings.Select(g => Math.Abs(g.X - nearby)).ToArray();
                int nearIndex = Array.IndexOf( nears,nears.Min());
                Point32 nearPeak = gratings[nearIndex];
                int pointsIndex = points.IndexOf(nearPeak.X);
                Point32[] subPoints = points.SubArrayPeriod(pointsIndex - 5, pointsIndex + 5);//取主瓣左右各5个点来计算峰值
                return FindMaxCurveFit(subPoints, xStep);
            }
            else
            {
                return FindMaxCurveFit(points, xStep);
            }
        }

        /// <summary>
        /// 用多项式拟合方法找阵列的最大值
        /// 注1：找出来的最大值是拟合出来的，通常与数列的最大值不相同
        /// 注2：提供的曲线最好是单峰的。
        /// </summary>
        /// <param name="points">要搜索的数组</param>
        /// <param name="xStep">x搜索步长</param>
        /// <returns>曲线拟合和的最大值</returns>
        public static Point32 FindMaxCurveFit(Point32[] points,double xStep=0.01)
        {//算法：1. 找到最大值 2. 寻找最大值左右两个点，共5个点用于曲线拟合 3.最大最小X值划分100份，求中间最大值

            var count = points.Count(pp => float.IsNaN(pp.Y));
            if(count>0)
            {
                throw new Exception("输入数据含有NaN值，无法参与计算！");
            }
            var maxPoint = FindMax(points);
            int indexMax = Array.IndexOf(points, maxPoint);
            if (indexMax == 0 || indexMax == points.Length - 1) return maxPoint;
            
            float[] xArray;
            float[] yArray;
            if (indexMax>=2 && indexMax < points.Length - 2)
            {
                yArray = new float[] { points[indexMax - 2].Y, points[indexMax - 1].Y, points[indexMax].Y, points[indexMax + 1].Y, points[indexMax + 2].Y };
                xArray = new float[] { points[indexMax - 2].X, points[indexMax - 1].X, points[indexMax].X, points[indexMax + 1].X, points[indexMax + 2].X };
            }
            else 
            {
                yArray = new float[] {  points[indexMax - 1].Y, points[indexMax].Y, points[indexMax + 1].Y };
                xArray = new float[] {  points[indexMax - 1].X, points[indexMax].X, points[indexMax + 1].X };
            }           
            double[] p = Fit.Polynomial(xArray.ToDouble(), yArray.ToDouble(), xArray.Length - 1);
            //多项式拟合的系数
            
            double deltaX = (xArray.Max() - xArray.Min()) / (100f );//步长缩短100倍去搜索最大值
            deltaX = deltaX < xStep ? deltaX : xStep;
            int N = (int)((xArray.Max() - xArray.Min()) / deltaX);
            double[] y = new double[N];
            for (int i = 0; i < N; i++)
            {
                double x = xArray.Min() + i * deltaX;
                y[i] = PolyValue(p, x);
            }
            double ymax = y.Max();
            double xMax = (float)(xArray.Min() + Array.IndexOf(y, ymax) * deltaX);
            return new Point32((float)xMax, (float)ymax);
        }
        /// <summary>
        /// 求解以p[]为系数的多项式，在x处的值
        /// </summary>
        /// <param name="p"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double PolyValue(double []p,double x)
        {
            double f = 0;

            for (int j = 0; j < p.Length; j++)
            {
                f += p[j] * Math.Pow(x, j);
            }
            return f;
        }
        /// <summary>
        /// 查找最大峰值
        /// 如果没有峰值，默认返回0，索引值返回-1
        /// </summary>
        /// <param name="array"></param>
        /// <param name="peakIndex">如果无峰值，返回-1</param>
        /// <returns>如果无峰值，返回0</returns>
        public static double FindMaxPeak(double[] array, out int peakIndex, float waveletFilter = 0)
        {
            double max = 0;
            double[] peaks = FindPeaks(array, out int[] peakIndexs, waveletFilter);
            if (peaks == null) peakIndex = -1;
            else
            {
                max = peaks.Max();
                int ind = IndexOf(peaks, max);
                peakIndex = peakIndexs[ind];
            }

            return max;
        }
        /// <summary>
        /// 寻找最大峰值
        /// 比较左右各order个值，均小于考察点视为峰值。
        /// </summary>
        /// <param name="points"></param>
        /// <param name="order">考察点的左右各取几个点来判断峰值，推荐2</param>
        /// <returns>如果成功返回峰值点（坐标和值）,否则返回PointF(float.MinValue, float.MinValue)</returns>
        public static Point32 FindMaxPeak(Point32[] points, int order=2)
        {
              var count = points.Count(pp => float.IsNaN(pp.Y));
            if (count > 0)
            {
                throw new Exception("输入数据含有NaN值，无法参与计算！");
            }
            List<Point32> peaks = FindPeaks(points, order);
            if (peaks.Count() > 0)
            {
                float max = peaks.Max(d => d.Y);
                return peaks.FirstOrDefault(d => d.Y == max);
            }
            else
            {
                return new Point32(float.MinValue, float.MinValue);
            }

        }
        /// <summary>
        /// 寻找数列的所有峰值点,适用于完整的周期性数据，查找峰值时会考虑边缘点是否为峰值
        /// 如某个角度范围内的方向图
        /// </summary>
        /// <param name="points"></param>
        /// <param name="order">阶数：两边各取几个点来参与判断</param>
        /// <param name="ignoreLevel">忽略与临近的峰值点相差小于此值的峰值</param>
        /// <returns></returns>
        public static List<Point32> FindPeaksPeriod(Point32[] points, int order = 2, double ignoreLevel = 0)
        {
              var count = points.Count(pp => float.IsNaN(pp.Y));
            if (count > 0)
            {
                throw new Exception("输入数据含有NaN值，无法参与计算！");
            }
            //var pp = SmoothPoints(points, 5);//做一次平滑处理，避免找错峰值点
            List<Point32> list = new List<Point32>();
            //for (int i = order; i < points.Length - order; i++)
            for (int i = 0; i < points.Length; i++)
            {

                //var subarr = MyMath.SubArray(points, i - order, i + order + 1);
                var subarr = points.SubArrayPeriod(i - order, i + order + 1);//
                var res = subarr.Where(a => a.X != points[i].X && a.Y >= points[i].Y);
                if (res.Count() == 0)//左右两边的y均小于中心值
                {
                    list.Add(points[i]);
                }
            }
            if (ignoreLevel > 0)//如果要求忽略掉某个电平内的副瓣
            {
                List<Point32> peaksToRemove = new List<Point32>();//等待删除的峰值
                for (int i = 1; i < list.Count; i++)
                {
                    if (Math.Abs(list[i].Y - list[i - 1].Y) < ignoreLevel)
                    {
                        Point32 toRemove = list[i].Y < list[i - 1].Y ? list[i] : list[i - 1];
                        peaksToRemove.Add(toRemove);
                    }
                }
                if (peaksToRemove.Count > 0)
                {
                    foreach (var pk in peaksToRemove)
                    {
                        list.Remove(pk);
                    }
                }
            }
            return list;
        }
        /// <summary>
        /// 寻找数列的所有峰值点。查找峰值点时，会忽略边缘order个点
        /// 适用于非周期性数据，如某个角度范围内的方向图
        /// </summary>
        /// <param name="points"></param>
        /// <param name="order">阶数：两边各取几个点来参与判断</param>
        /// <param name="ignoreLevel">忽略与临近的峰值点相差小于此值的峰值</param>
        /// <returns></returns>
        public static List<Point32> FindPeaks(Point32[] points, int order = 2, double ignoreLevel=0)
        {
              var count = points.Count(pp => float.IsNaN(pp.Y));
            if (count > 0)
            {
                throw new Exception("输入数据含有NaN值，无法参与计算！");
            }
            //var pp = SmoothPoints(points, 5);//做一次平滑处理，避免找错峰值点
            List<Point32> list = new List<Point32>();
            for (int i = order; i < points.Length - order; i++)
            {

                var subarr = MyMath.SubArray(points, i - order, i + order + 1);
                //var subarr = points.SubArrayPeriod(i - order, i + order + 1);//
                var res = subarr.Where(a => a.X != points[i].X && a.Y >= points[i].Y);
                if (res.Count() == 0)//左右两边的y均小于中心值
                {
                    list.Add(points[i]);
                }
            }
            if (ignoreLevel > 0)//如果要求忽略掉某个电平内的副瓣
            {
                List<Point32> peaksToRemove = new List<Point32>();//等待删除的峰值
                for (int i = 1 ; i < list.Count; i++)
                {
                    if(Math.Abs( list[i].Y-list[i-1].Y)<ignoreLevel)
                    {
                        Point32 toRemove = list[i].Y < list[i - 1].Y ? list[i] : list[i - 1];
                        peaksToRemove.Add(toRemove);
                    }
                }
                if(peaksToRemove.Count>0)
                {
                    foreach (var pk in peaksToRemove)
                    {
                        list.Remove(pk);
                    }
                }
            }
            return list;
        }
        /// <summary>
        /// 寻找数列的所有峰值点
        /// 将数列视为首尾相接搜索
        /// </summary>
        /// <param name="points"></param>
        /// <param name="order">阶数：两边各取几个点来参与判断</param>
        /// <returns></returns>
        public static List<Point32>SearchPeaksPeriod(Point32[] points,int order=2)
        {
              var count = points.Count(pp => float.IsNaN(pp.Y));
            if (count > 0)
            {
                throw new Exception("输入数据含有NaN值，无法参与计算！");
            }
            List<Point32> list = new List<Point32>();
            for (int i = 0; i < points.Length; i++)
            {
                var subarr = points.SubArrayPeriod(i - order, i + order+1);
                var res = subarr.Where(a => a.X != points[i].X && a.Y >= points[i].Y);
                if(res.Count()==0)//左右两边的y均小于中心值
                {
                    list.Add(points[i]);
                }
            }
            
            return list;
        }
        /// <summary>
        /// 用曲线拟合方式寻找峰值
        /// </summary>
        /// <param name="points"></param>
        /// <param name="curveFit">是否采用曲线拟合方式</param>
        /// <returns></returns>
        public static List<Point32> FindPeaks(Point32[] points,  bool curveFit)
        {
              var count = points.Count(pp => float.IsNaN(pp.Y));
            if (count > 0)
            {
                throw new Exception("输入数据含有NaN值，无法参与计算！");
            }
            if (points.Length < 3)
            {
                throw (new Exception("没有足够数据(至少3points)用以寻找峰值"));
            }
            List<Point32> peaks = new List<Point32>();
            for (int i = 1; i < points.Length - 1; i++)
            {
                if (points[i].Y > points[i - 1].Y  && points[i].Y > points[i - 1].Y )
                {
                    if (curveFit)
                    {
                        Point32[] subPoints;
                        if (points.Length >= 5 && i >= 2 && i <= points.Length - 3)
                        {//做5项式拟合
                            subPoints = SubArray(points, i - 2, i + 2);
                        }
                        else
                        {//3项式拟合
                            subPoints = SubArray(points, i - 1, i + 1);
                        }
                        Point32 peak = FindMaxCurveFit(subPoints);
                        peaks.Add(peak);
                    }
                    else
                    {
                        peaks.Add(points[i]);
                    }
                }
            }
            return peaks;
        }
        /// <summary>
        /// 获取子数列。
        /// 如果lowIndex小于 upIndex,则取其内区间子数列
        /// 否则，取外区间子数列（原数列中array[upIndex...Length-1,0... lowIndex]
        /// </summary>
        /// <typeparam name="T"></typeparam>
        /// <param name="array"></param>
        /// <param name="lowIndex">下边界</param>
        /// <param name="upIndex">上边界</param>
        /// <returns></returns>
        public static T[] SubArray<T>(T[] array, int lowIndex, int upIndex)
        {
            List<T> sa = new List<T>();
            if(lowIndex<=upIndex)
            {
                for (int i = lowIndex; i < upIndex; i++)
                {
                    sa.Add(array[i]);
                }
            }
            else
            {
                for (int i = lowIndex; i < array.Length; i++)
                {
                    sa.Add(array[i]);
                }
                for (int i = 0; i < upIndex; i++)
                {
                    sa.Add(array[i]);
                }
            }
            
            return sa.ToArray();
        }
       
        public static T[,] SubArray<T> (T[,] array,int lowIndex,int upIndex)
        {
            int row = upIndex - lowIndex + 1;
            int col = array.GetLength(1);
            T[,] sub = new T[row, col];
            for (int i = 0; i < row; i++)
            {
                for (int j = 0; j < col; j++)
                {
                    sub[i, j] = array[i + lowIndex, j];
                }
            }
            return sub;
        }
        /// <summary>
        /// 将pointf阵列中x的元素移动到中心位置
        /// </summary>
        /// <param name="points"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        public static Point32[] MovePointsToMidle(Point32[] points,double x)
        {
            int index = points.IndexOf(x);
            return MoveToMiddle(points, index);

        }
       
        /// <summary>
        /// 将阵列元素位移到以指定索引为中心
        /// </summary>
        /// <typeparam name="T"></typeparam>
        /// <param name="array"></param>
        /// <param name="midIndex"></param>
        /// <returns></returns>
        public static T[] MoveToMiddle<T>(T[]array,int midIndex)
        {
            int N = array.Length;
            int mid = N / 2;
            int m = midIndex;
            
            T[] arr = new T[array.Length];
            if(m<mid)
            {
                for (int i = 0; i <N; i++)
                {
                    if (i < mid - m)
                        arr[i] = array[i + mid + m ];
                    else 
                        arr[i] = array[i + m - mid];
                }
            }
            else
            {
                for (int i = 0; i < N; i++)
                {
                    if (i < N+mid-m)
                        arr[i] = array[i + m - mid];
                    else 
                        arr[i] = array[i + m - mid - N];
                   
                }
            }

            return arr;
        }
        /// <summary>
        /// 将数列按指定的方向旋转rotate量         
        /// 本方法假定数列是首位相接的。旋转时，数列的索引将按指定方向做位移
        /// </summary>
        /// <typeparam name="T"></typeparam>
        /// <param name="dataArray">要旋转的数列</param>
        /// <param name="rotate">旋转的索引数</param>
        /// <param name="antiClockwise">是否逆时针旋转</param>
        /// <returns></returns>
        public static T[] RotateArray<T>(T[] dataArray, int rotate, bool antiClockwise = true)
        {
            if(dataArray.Length<rotate)
            {
                rotate -= dataArray.Length;
            }
            T[] array = new T[dataArray.Length];
            if (antiClockwise)
            {
                for (int i = 0; i < array.Length - rotate; i++)
                {
                    array[i] = dataArray[rotate + i];
                }
                for (int i = array.Length - rotate; i < dataArray.Length; i++)
                {
                    array[i] = dataArray[i - (array.Length - rotate)];
                }
            }
            else
            {
                for (int i = 0; i <= rotate; i++)
                {
                    array[i] = dataArray[array.Length - rotate + i - 1];
                }
                for (int i = rotate + 1; i < dataArray.Length; i++)
                {
                    array[i] = dataArray[i - rotate - 1];
                }
            }
            return array;
        }

        /// <summary>
        /// 从二维矩阵中提取指定列
        /// </summary>
        /// <typeparam name="T"></typeparam>
        /// <param name="array">array[,]</param>
        /// <param name="colIndex"></param>
        /// <returns></returns>
        public static T[] SpecifyColsFrom2DArray<T>(T[,] array,int colIndex)
        {
            T[] data = new T[array.GetLength(0)];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = array[i, colIndex];
            }
            return data;
        }
        public static double[] RealPart(Complex[] complex)
        {
            double[] r = new double[complex.Length];
            for (int i = 0; i < r.Length; i++)
            {
                r[i] = complex[i].Real;
            }
            return r;
        }
        public static double[] ImagPart(Complex[] complex)
        {
            double[] r = new double[complex.Length];
            for (int i = 0; i < r.Length; i++)
            {
                r[i] = complex[i].Imaginary;
            }
            return r;
        }
        public static double[] Abs(Complex[] complex)
        {
            double[] r = new double[complex.Length];
            for (int i = 0; i < r.Length; i++)
            {
                r[i] = complex[i].Magnitude;
            }
            return r;
        }
        public static double[] Phase(Complex[] complex)
        {
            double[] r = new double[complex.Length];
            for (int i = 0; i < r.Length; i++)
            {
                r[i] = complex[i].Phase;
            }
            return r;
        }
        public static double[] Vswr(Complex[] S11)
        {
            double[] vs = new double[S11.Length];
            for (int i = 0; i < vs.Length; i++)
            {
                vs[i] = (1 +  S11[i].Magnitude) / (1 - S11[i].Magnitude);
            }
            return vs;
        }

        

        public static double[] Db(Complex [] sp)
        {
            double[] db = new double[sp.Length];
            for (int i = 0; i < sp.Length; i++)
            {
                db[i] = 20.0 * Math.Log10(sp[i].Magnitude);
            }
            return db;
        }
       
        /// <summary>
        /// 寻找指定数据点的最大值
        /// </summary>
        /// <param name="subPats"></param>
        /// <returns>用PointF型返回最大值点坐标及峰值</returns>
        public static Point32 FindMax(Point32[] points)
        {
            var maxy = points.Max(p => p.Y);
            var x = points.FirstOrDefault(p => p.Y == maxy).X;
            return new Point32(x, maxy);
            
        }
        /// <summary>
        /// 寻找指定数据点的最大值
        /// </summary>
        /// <param name="points"></param>
        /// <returns></returns>
        public static Point32 FindMin(Point32[] points)
        {
            var minY = points.Min(p => p.Y);
            var x = points.FirstOrDefault(p => p.Y == minY).X;
            return new Point32(x, minY);

        }

    }
}
